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Real-Time Simulation of Non-Equilibrium Transport of Magnetization 
in Large Open Quantum Spin Systems driven by Dissipation 
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Using quantum Monte Carlo, we study the non-equilibrium transport of magnetization in large 
open strongly correlated quantum spin \ systems driven by purely dissipative processes that conserve 
the uniform or staggered magnetization, disregarding unitary Hamiltonian dynamics. We prepare 
both a low-temperature Heisenberg ferromagnet and an antiferromagnet in two parts of the system 
that are initially isolated from each other. We then bring the two subsystems in contact and study 
their real-time dissipative dynamics for different geometries. The flow of the uniform or staggered 
magnetization from one part of the system to the other is described by a diffusion equation that 
can be derived analytically. 
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Simulating the real-time evolution of large strongly 
correlated quantum systems is notoriously difficult, due 
to the dimension of the Hilbert space, which grows ex¬ 
ponentially with the system size. In this case, Monte 
Carlo methods are usually not applicable because impor¬ 
tance sampling is prevented by severe sign or complex 
phase problems [1]. While in Euclidean time some se¬ 
vere sign problems have been solved using the meron- 
cluster algorithm Hi or the fermion bag method [3-EJ , 
until recently real-time simulations of quantum systems 
have been limited to small volumes that are accessible 
to exact diagonalization, or to gapped 1-d systems to 
which the time-dependent density matrix renormaliza¬ 
tion group (DMRG) [3,0 can be applied. Even then, due 
to the growth of entanglement, only moderate time inter¬ 
vals can be investigated UGH- Dynamical phenomena in 
non-equilibrium quantum systems have been studied in 
Mz]. Recently, we have developed a new Monte Carlo 


method that allows us to simulate the real-time evolution 
of large strongly coupled quantum systems in any dimen¬ 
sion for an arbitrary amount of time, for specific dynam¬ 
ics driven by purely dissipative processes that are de¬ 
scribed by a Lindblad equation 28,[29]. In particular, the 
unitary time-evolution driven by a Hamiltonian, which 
would give rise to a severe complex phase problem, has 
been replaced by a dissipative process. Still severe sign 
problems arise even for the purely dissipative dynam¬ 
ics, but they have been solved analytically by identifying 
exact cancellations in the corresponding real-time path 
integral. Purely dissipative processes play an important 
role in quantum information processing, for example, in 
order to prepare specific states for quantum computation 
30l [ffij or entanglement generation [ 


The control of 


quantum systems by measurements has been investigated 
36, 37) • Ultracold atoms in optical lattices or trapped 


m 


ions provide platforms in which such dynamics can be 
engineered in quantum simulation experiments (38144^. 

In this paper, our primary goal is not yet to make con¬ 


tact with concrete cold atoms experiments. Instead, we 
demonstrate that our ability to classically simulate the 
real-time dynamics of engineered dissipative processes in 
large open quantum spin systems puts us in a unique po¬ 
sition to study transport phenomena far away from equi¬ 
librium. Such processes thus provide a bridge between 
classical and quantum simulations of real-time quantum 
dynamics. Here we investigate a low-temperature Heisen¬ 
berg ferromagnet and an antiferromagnet which are ini¬ 
tially isolated from each other in two separate parts of 
the volume. The two parts, which act as large reser¬ 
voirs of uniform or staggered magnetization, are then 
put in contact and evolve in time according to a dissi¬ 
pative process which either conserves the uniform or the 
staggered magnetization. The corresponding conserved 
quantity then flows from its reservoir into the other half 
of the system, through an opening whose size we vary. 
The non-equilibrium diffusive processes are driven by the 
gradient of the corresponding conserved quantity. They 
come to an end only when the staggered or uniform mag¬ 
netization is homogeneously distributed throughout the 
entire system. Remarkably, certain aspects of the dy¬ 
namics are described by a classical diffusion equation 
which can be derived analytically from the underlying 
dissipative quantum dynamics. Significantly extdending 
previous work [H, Q , the current setting allows us to 
study the diffusion process of the conserved quantity in 
real-space. 


We consider systems of quantum spins ^ on a square 
lattice, which are dissipatively coupled to their envi¬ 
ronment. The dynamics is characterized by a set of 
Lindblad operators L 0k , that obey (1 — 57)1 + 

o k Li k L 0k = 1, where 5 is a small time-step. The 
Lindblad operators induce quantum jumps and 7 deter¬ 
mines their probability per unit time. We will analyt¬ 
ically derive the relation between the parameter 7 and 
the diffusion coefficient of the classical diffusion equa¬ 
tion. The time-evolution of the density matrix is then 
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determined by the Lindblad equation 
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We will consider two different dissipative processes 
whose jump operators L 0k = yJe^P 0k are determined 
by operators Po k that project on the eigenstates of an 
observable O with eigenvalue o^. For the first process 
(process 1), which conserves the uniform magnetization 
vector, the observable is the total spin O ^ = (S x + S y ) 2 
of a pair of spins S x and S y located on neighboring lattice 
sites x and y. The projection operators corresponding to 
total spin 1 or 0 are then given by 
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, the conservation of the 
total spin in this dissipative process implies that the 
low-momentum modes of the magnetization equilibrate 
very slowly. The second dissipative process (process 2), 
which conserves the 3-component of the staggered mag¬ 
netization, is characterized by the observable O^ = 


S+S++S- 


S with the three projection operators 


(\ 
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In this case, as we showed in [29|, the high-momentum 
modes of the magnetization (namely those with momenta 
near the conserved (7r, 7r)-mode representing the stag¬ 
gered magnetization) equilibrate very slowly. Both dissi¬ 
pative processes ultimately converge to a trivial infinite- 
temperature density matrix that is proportional to the 
unit-matrix, at least within the sector defined by the 
value of the corresponding conserved quantity. 

As discussed in detail in 28|, [29], the Lindblad equa¬ 
tion can be represented by a path integral consisting of 
a Euclidean time contour that defines an initial density 
matrix in thermal equilibrium and a real-time Schwinger- 
Keldysh contour [45|, |46| that leads from an initial time 
to to a final time tf and back. Remarkably, the proba¬ 
bility to reach a specific final state |/) can be computed 
very efficiently with a loop-cluster alg orithm, similar to 
the one used in Euclidean time |47449] . The cluster rules 
have been discussed in detail in [29|. 

We consider a spin | Heisenberg model with Hamil¬ 


tonian H = J ^2 { 


xy) k 


‘Sy. 


In order to prepare an ini¬ 


tial density matrix we consider an L x 2L lattice that 


is divided into two subsystems of size L x L each with 
individual periodic boundary conditions. One system is 
antiferromagnetic (with J > 0) and the other is ferro¬ 
magnetic and has the opposite exchange coupling. Both 
subsystems are initialized at the same temperature. The 
initial density matrix is then subjected to one of the two 
dissipative real-time processes. During the real-time pro¬ 
cess the two subsystems are put in contact through two 
openings of size L' < L on opposite sides of both sys¬ 
tems. This is achieved by changing the original boundary 
conditions with period L on two sets of L' links. These 
links connect the two subsystems, so that the total sys¬ 
tem now has boundary conditions with period 2 L in a 
strip of transverse size L' and the original pair of bound¬ 
ary conditions with period L on the remaining strip of 
transverse size L — L '. The transverse direction of size L 
always has ordinary periodic boundary conditions. Using 
the loop-cluster algorithm we calculate the expectation 
value of the 3-component for each spin at the time t / 
when the 3-components of all spins are finally being mea¬ 
sured. The data are separately analyzed for each total 
value of the conserved uniform or staggered magnetiza¬ 
tion. By using an improved estimator similar to the one 
constructed in [50|, [5_l| , we increase the statistics by a 
factor that grows exponentially with the number of loop- 
clusters. This improves the accuracy of the numerical 
data very substantially and leads to the results depicted 
in Fig. [T] (uniform magnetization) and Fig. [2] (staggered 
magnetization). As we have discussed in detail in (28. 29j|. 
the dissipative processes give rise to different time scales. 
While process 1 quickly destroys the initial antiferromag¬ 
netic order over a time scale l/y, the conserved uniform 
magnetization undergoes a much slower diffusion pro¬ 
cess. In particular, in process 1 the magnetization modes 
with low momentum p equilibrate only over time scales 
l/(ya 2 p 2 ), where a is the lattice spacing. Similarly, in 
process 2, which conserves the staggered magnetization, 
the modes with momenta near ( 7 r, 7 r) are severely slowed 
down. The dissipative dynamics can be characterized as 
a heating process that affects different modes at different 
time scales. While the underlying diffusive processes 
are quantum mechanical, the resulting expectation val¬ 
ues of the conserved uniform or staggered magnetization 
are described by a classical diffusion equation 

dtpx(t) = | E I Px+aM ~ 2 Px(t) + P x -ai(t)] ■ (4) 


Here p x ( t ) is the expectation value of the conserved quan¬ 
tity at the lattice site x at time t and i is the unit-vector 
in the i-direction. Interestingly, the classical diffusion 
equation can be derived analytically from the underly¬ 
ing quantum spin dynamics, and the diffusion coefficient 
is determined by the parameter 7 that drives the Lind¬ 
blad process of eq. ©• The lattice diffusion equation (g]) 
results from the continuity equation 
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FIG. 1: [Color online] Real-time evolution of the uniform 
magnetization on a 32 x 64 lattice with an opening of size 
L' = 4a for a total uniform magnetization value M u — 
\(L/a) 2 = 512 at initial temperature /3J = 80. Typical con¬ 
figurations (left) and expectation values of the uniform mag¬ 
netization (right) at time t — 0 (top), 50/7 (middle), and 
500/7 (bottom). 


FIG. 2: [Color online] Real-time evolution of the staggered 
magnetization on a 32 x 64 lattice with an opening of size 
L' = 4a for a total staggered magnetization value M s = 
| (L/a) 2 = 384 at initial temperature /3J — 80. Typical con¬ 
figurations (left) and expectation values of the staggered mag¬ 
netization (right) at time t — 0 (top), 50/7 (middle), and 
500/7 (bottom). 
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FIG. 4: [Color online] Real-time evolution of the total uniform 
magnetization in the first subsystem M for different values of 
L' (L = 32 a, total uniform magnetization M u = L/a ) 2 — 
512/. Inset: Late-time relaxation rate as a function of L'. 


FIG. 3: Configurations of two neighboring spins evolving in 
time, together with the resulting values for p x (t) and j x ,i(t) for 
the dissipative process 1 and 2, that conserves the uniform and 
staggered magnetization, respectively. The current is driven 
by the gradient of the corresponding density. 


dtpx(t) + ^ ix-alp^) 


= 0, 


(5) 


combined with the lattice gradient equation 

= —Y [p x+a l(t) - Px{t)} . (6) 

Here ) is the conserved (uniform or staggered) mag¬ 
netization current density that flows from the lattice site 
x to the neighboring lattice site x + ai at the time t. The 
continuity equation (|5j) and the gradient equation (|6j) can 
be derived from the underlying real-time path integral 
that was discussed in detail in [28|, [29]. The correspond¬ 
ing spin configurations together with the resulting values 
for p x (t) and j x ,i(t ) are illustrated in Fig. [3] for the two 
dissipative processes. 

We have also investigated the time-dependence of the 
total uniform magnetization in the first subsystem (ini¬ 
tially ferromagnetic) as a function of the opening size L' 
in dissipative process 1 (cf. Fig. |4|). The final state, for 
which the magnetization is homogeneously distributed 
throughout the entire system, is reached exponentially 
at long times. The relaxation rate then depends linearly 
on the opening size L' over a wide range of values of L'. 

For the largest possible size of openings, L' = L, the 
diffusion equation reduces to a 1-d problem which can 
even be solved analytically. The resulting profile of the 
magnetization density, illustrated in Fig. [5j is given by 


p (t) = Po lasin (ff (2^ +a)) 

2 tzi 

(n odd) 


Lsinfff) 


x exp (-2 tW (=£) t) + f . (7) 


Certain features of the dissipative processes discussed 
here resemble classical physics. For example, if finally 
all spins are projected along the 3-axis at the end of the 
real-time evolution, the spin configurations on the two 
branches of the Keldysh contour become identical 28|, 29] 
and their evolution reduces to a Kawasaki dynamics [52] 
(cf. Fig. [3]), which can be captured by a classical diffusion 
equation. Other aspects of the same dynamics, includ¬ 
ing the time-evolution of entanglement, do not have this 
feature, thus underscoring the quantum nature of the 
corresponding real-time processes. In particular, while 
the expectation value of the conserved quantity obeys a 
classical diffusion equation, its probability distribution 
can only be calculated quantum mechanically. We em¬ 
phasize that the presented method is not restricted to 
probing only diagonal elements of the density matrix. 
Most notably, two-point correlation functions reflecting 
off-diagonal entries of the density matrix could also be 
measured very efficiently via improved estimators [53| . 



FIG. 5: [Color online] The 1-d profile of the uniform magne¬ 
tization density (L' = L = 32 a) evolves from a step function 
at the initial time to a uniform distribution at late times. 
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It would be most interesting to investigate the non- 
dissipative pure Hamiltonian dynamics of large closed 
quantum systems. Due to very severe complex phase 
problems this is most likely impossible on a classical com¬ 
puter. On the other hand, quantum simulators, for exam¬ 
ple, using ultracold atoms in optical lattices, are ideally 
suited for such investigations. It is conceivable to exper¬ 
imentally design a dissipative environment which acts as 
a projector on singlet and triplet states (process 1) with 
current technology Q, whereas the realization of process 
2 is probably more involved. On the other hand, as we 
have shown, engineered purely dissipative processes are 
accessible to very efficient real-time simulation of large 
open quantum systems using classical computers. Such 
real-time processes thus provide a bridge between clas¬ 
sical and quantum simulation. It will be most interest¬ 
ing to explore other processes, including a weakly cou¬ 
pled Hamiltonian or non-Hermitean Lindblad operators, 
in order to explore the territory connecting classical and 
quantum simulation of quantum dynamics in real time. 
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